# price update formula for fixed point iteration
price.update <- function(price,nu) {
  #market shares
  pim <- prob.mat(price)
  eps <- -ownelas(pim)
  r <- cannib(pim,price)
  # apply the first order condition for profit-max.
  p.star <- mc * (eps+1)/(eps-r)
  # weight the nu price
  nu*p.star +(1-nu)*price
}
#
#prob.mat creates the whole prob matrix 
prob.mat <- function(p) {
  # matrix of indirect utilties: row= household, col = variety
  ubar <- b.h %*% t(x) - alpha.h %*% t(log(p)) +rep(1,n.hh) %*% t(xi)
  # MATRIX [n.hh X n.models] of household-level probabilities of variety j
#  exp(ubar)/(U0+rowSums(exp(ubar)))  # U0= 0 or 1 in global env.
  ubar <- cbind(ubar,rep(log(U0),n.hh))
  # MATRIX [n.hh X n.models] of household-level probabilities of variety j
  I <- matrixStats::rowLogSumExps(ubar)
  return(exp(ubar-I)[,1:n.models])
  }
#  
#elasticities of OWN share wrt OWN price (an n.models length vector)
ownelas <- function(prm) {
  colSums(-alpha.h*(1-prm)*omega.mat(prm)) # summing over individual weighted elasticties
}
#  
#matrix of omega for own elasticity calculation (matrix of n.hh X n.models).
omega.mat <- function(prm) {
  omega.mat <- matrix(0,nrow=n.hh,ncol=n.models)
  for(i in 1:n.models) {
    omega.mat[,i] =omega(prm,i) 
  }
  return(omega.mat)
}
#  
#derivatives of OTHER CO-OWNED shares wrt OWN price (an n.models length vector)
cannib <- function(prm,price)  {
  s <- share.mces(prm)
  S <- s*Y
  D.cross <- crossDeltaMCES(MM=co.own,HH=prm,A=alpha.h,S=S,y=y.h)
#  D.cross <- build.DeltaX(prm)
  L <- (price-mc)/price
   return((D.cross %*% (L*s))/s )  # called "r" in theory
}
#  
#matrix of cross price elasticities, changing price of m (row), impact on j(cols)
build.DeltaX <- function(prm) { # in prm rows are hh, co.own mXm plucked from env.
  Delta <- matrix(0,nrow=n.models,ncol=n.models)
  for(m in 1:n.models) {
    for(j in 1:n.models) {
      if(co.own[m,j]) Delta[m,j] = sum(alpha.h*prm[,m]*omega(prm,j))  # epsilon_{jm} * Omega_{jm} in theory
    }
  }
  return(Delta)
}
#
# n.hh length vector of sales-shares for model m. (nothing to do with ownership matrix)
omega <- function(prm,m) {
  (prm[,m]*y.h)/sum(prm[,m]*y.h)  
}
#  
#variety-level market shares 
share.mces <- function(prm) {
  colSums(prm*y.h)/Y  
}
#variety-level average income of consumers   
avg.inc <- function(prm) {
  num <- colMeans(prm*y.h) 
  denom <- colMeans(prm)
  num/denom
}  
